function RR=mconv(M,N,m,n);

% MCONV    subroutine to produce the product of two matrix polynomials.
%          R(L) = M(L) * N(L)
%             where M(L) = M0 + M1 L^1 + .... Mm L^m
%                   N(L) = N0 + N1 L^1 + .... Nn L^n
%          the polyomial coefficients are stored in the input matrices
%             M=[M0 M1 .... Mm] and N=[N0 N1 .... Nn]
%          and the output matrix
%             RR=[R0 R1 ... Rr]  with r=n+m.

% The strategy in this program is to use multiplication of ever expanding
% matrices to produce the convolution coefficients.
%   i = 0  RR0 = M0*N0
%   i = 1  RR1 = MC*NC = [M0 M1]* [N1  = M0*N0 + M1*N1
%                                  N0]
%   etc.  For values of i exceeding m(n) the MC(NC) matrix is padded with zeros.
% This strategy is based on the fact that matrix multiplication is fast in 
% most computer programs while loops are slow.

dM=size(M);

rm=dM(1);        % number of rows of Mi matrix;
cm=dM(2)/(m+1);  % number of columns of Mi matrix;

dN=size(N);
rn=dN(1);        % number of rows of Ni matrix;
cn=dN(2)/(n+1);  % number of columns of Ni matrix;

if (cm~=rn)
  disp('matrices illspecified: mconv terminating')
  RR=0;
else
  % initialize output matrix.
  RR=zeros(rm,cn*(m+n+1));
  % Create matrices of coefficients padded with zeros. 
  MM=zeros(rm,cm*(m+n+1));
  MM(:,1:dM(2))=M;
  NN=zeros(rn,cn*(m+n+1));
  NN(:,1:dN(2))=N;
  % Initialize current matrices
  MC=MM(:,1:cm);
  NC=NN(:,1:cn);
  i=0;
  while i<= (n+m-1);
   RR(:,i*cn+1:i*cn+cn)=MC*NC;
   MC=[MC MM(:,(i+1)*cm+1:(i+2)*cm)];  % update MC matrix
   NC=[NN(:,(i+1)*cn+1:(i+2)*cn)
       NC];     % update NC matrix
   i=i+1;
  end
  RR(:,i*cn+1:i*cn+cn)=MC*NC;
end

